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Abstract 



The thermal lattice BGK model is a recently suggested numerical tool which aims 
to simulate thermohydrodynamic problems. We investigate the quality of the lattice 
BGK simulation by calculating temperature profiles in the Couette fiow under difli"erent 
Eckert and Mach numbers. A revised lower order model is proposed and the higher 
order model is proved once again to be advantageous. 

1 Introduction 

The lattice BGK method can be viewed as the latest development of the lattice Boltzmann(LB) 
method[2], which is a derivation of the lattice gas automata(LGA) model[l] for the simulation of 
fiuid dynamics. Such a development was achieved by introducing the elaborate Bhatnagar-Gross- 
Krook collision operator[6] into the lattice Boltzmann equation(LBE)[3][4][5]. Although this single 
time relaxation approximation(STRA) made on the collision term of the discrete kinetic equation 
looks like oversimplified, the lattice BGK models amazingly reproduced well the complexities of 
fiuid fiows[7][8][9]. 

Recently, lattice BGK models in which thermal effects are included were explored by several 
authors[10][ll][13]. These models are established by treating the conservation of particle kinetic 
energy non-trivially, which can be realized with the use of multi-speed particle distributions. 
The employment of the composed lattice is basically required, e.g. Alexander[10] employed the 
composed hexagonal lattice in two dimensional space, while others used the composed square 
lattice in the D dimensional space. One may further divide these thermal lattice BGK models 
into two classes: the lower order models(Alexander and Qian's models[ll]) and the higher order 
model(Chen's model[13]). The word "order" here refers to the order of the fiow speed u in its 
expansion of the equilibrium particle distribution, which appears in the lattice BGK equation as. 



Here, Npj^i represents the particle distribution on the ith link of the pk sub-lattice, Cpj^i is the link 
vector and therefore the vector of the particle fiight velocity, and r is the relaxation time period 
during which particle distributions would approach to their equilibrium values. The equilibrium 
particle distribution, A'^'^f , is written in the form of low speed expansion up to the 2nd(or 3rd) 
order of u for the lower order models and 4th order for the higher order model. Correspondingly, 
different levels of symmetries of the underlying lattice are required for these models, so that 
the momentum and heat fiux of the modeled fiuid can be expressed in an isotropic form in the 
macroscopic limit. Specifically, the nth rank particle velocity-moment tensor defined as 



is required to have an isotropic form up to n = 4 for the lower models and n = 6 for the higher order 
model. The latter requirement would prevent using the composed hexagonal lattice in 2D space 
and the FCHC lattice in the AD space, for that the highest rank of the isotropic velocity-moment 
tensor is four for both cases[14]. 

Needless to say the higher order thermal lattice BGK model is more accurate, which was proved 
by the numerical measurement of the decaying rates of fiow kinetic energy under different Mach 
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numbers[13]. In that case, the higher order model was shown to be free from the deviation effects 
caused by the nonlinear terms hidden in the r.h.s. of the momentum equation of the modeled 
fluid. In this study, however, we shall concentrate on the aspect of heat transfer modeling for the 
thermal lattice BGK models. The investigation will be carried out by measuring and comparing 
the temperature profiles in the Couette fiows under different Eckert numbers and Mach numbers. 
Models based on the composed square lattice are used in all these calculations. The key issues 
of the thermal lattice BGK model will be briefed in the next section, and numerical experiments 
and their results will be described subsequently. The lower and higher order models used in this 
paper will be labeled as follows, 

• 2D13VQ: Qian's lower order model, based on the 2D 13-link composed square lattice. 

• 2D13VC: a revised lower order model, based on the 2D 13-link composed square lattice. 

• 2D16V : Chen's higher order model, based on the 2D 16-link composed square lattice. 

Hence the discussions are confined in two dimensional cases, though the extension to three dimen- 
sional cases is straightforward. Some concluding remarks will be given in the last section. 

2 Thermal lattice BGK models 

2.1 Coordinates and symmetries of the composed lattice 

The coordinates of the base vectors of the square sub-lattices may be written as, 

D components in D dimensions 
k( ±1,±1,...,±1,0,...,0,0 ), (3) 
p components 

and its permutations. The number of non-zero components is p so that the moduli of such a vector 
is \cpki \ = k^/p. It is important to know, in the process of hydrodynamic derivation, the expression 

of the velocity-moment tensor T^^^ ^ which consists of the n-product of base vectors. The odd 
rank tensors vanish naturally by the definition itself. The even rank tensors, specifically 2nd, 4th 
and 6th rank tensors, can be written as follows in the D dimensional space, 

T!,lip = ^Pk^cp, (4) 

Tpkl/ljSCe = ^pk^al3~f6C^ + ^pk{6al3^~f6C^ + ...)+ ^pki^apT^t^SCe + •••)• 

Here, 6 is the Kronecker tensor and T is its higher order version. The ellipsis ". . ." stands for 
terms which can be obtained by permuting the indices of the foregoing term. The numerical 
values of parameters, such as i9pk, ippk, ^pk are listed in references[14][12] for one, two and 
three dimensions. The way to make the macroscopic fiux tensors isotropic is to tune the particle 
populations on different sub-lattices properly so that the sum of anisotropic parts, parts that are 
related to tppk, ^pk and f^p^, vanishes as a total effect. 



2.2 Equilibrium particle distribution 

When the fiow speed of the modeled fiuid is controlled to be much smaller than the fiight speed of 
particles, the local equilibrium particle distribution can be expanded around the uniform equilib- 
rium state. Considering fiexibilities in the residence and the number of particles, and the parity 
invariance of the square lattice, the low speed expansion is written as follows, 

^^ki = \k + Mpj,{CpMaUa) + GpkU^ + JpkiCpkiaUaf + (5) 
Qpk(y^pkia'^a)'^ H" ^pk(y^pkia'^a) H" ^pk(y^pkia'^a) ^ H" ^pk"^ H" ^(^U 'j . 
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For the lower order models, terms whose orders are higher than may be cut out, or one of 
the third order terms can be retained to make the energy equation more accurate on the Euler 
level[12]. Parameters of the expansion depend on the local conserved quantities, namely density 
p and thermal energy e. This dependence can be written in the following form, 

2 

Xpk = p'^Xpj.ie' . (6) 
;=o 

Here, Xpk may represent any one of Apk, Mpk, Spk so that Xpki is actually Gpki, rupuu Spki- 
Various constraints, which may ensure the definition of the conserved quantities, the vanishing 
of anisotropic parts in the macroscopic fiux tensors, and the nonexistence of unphysical artifacts 
in the macro-dynamic equations, can be imposed on these parameters. The number of such 
constraints is usually smaller than that of the parameters, so that the specification of a thermal 
lattice BGK model always involves some arbitrariness[ll][12]. This implies that some optional 
constraints may be employed either to minimize the anisotropic effects or to improve the accuracy 
of the model, which would be made clear in the following sections. 



2.3 Macroscopic flux tensors and transport coeflRcients 

In the hydrodynamic derivation for the lattice BGK models, the discrete kinetic equation shall 
be first Taylor expanded into a continuous form in the long wavelength and low frequency limits. 
The macro-dynamic equations for the conserved quantities can be obtained subsequently by using 
the multi-scale technique, which is a perturbative formulation method. The small quantity for 
perturbation is proportional to the local Knudsen number and is denoted as e. Correct forms of 
the macro-dynamic equations are ensured by correct expressions of the macroscopic fiux tensors, 
which are again ensured by the aforementioned constraints imposed on parameters of the low 
speed expansion of the equilibrium particle distribution. The resulted momentum fiux tensors on 
the Euler(e) order and Navier-Stokes(e^) order are given as, 

^«/3 = '^^p°k]cpkiaCpkil3 = J^PeSajS + pUaUjj , (7) 
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Here A'^^] is the equilibrium part and A'^^,- is the nonequilibrium part of the particle distribution. 
Heat fiux vectors may be expressed on different orders as follows. 
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Both the transport coefficients and the state equation of the modeled fiuid can be easily identified 
from these formulas, e.g. the shear viscosity and the heat conductivity read respectively as, 

2,1, 
P = DP^y'^~2'' 

(D + 2) ,1, 
= -^P<^-2^- (12) 
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Note that the accuracies of Eqs. (7), (8), (9) and (10) are different for models of different orders. 
For the lower order models, hJ^"^ and are accurate up to eu^ and et/^ orders, and hJ^"'^ and 

qa^ up to e^u^ and e^u orders. For the higher order model, the accuracies of the corresponding 
quantities are upgraded to eu'' , eu^ , e^u^ and e^u^ orders. 

3 Numerical investigations 

3.1 Couette flow 

The flow system is an extremely simple one, consisting of one moving boundary, one rest boundary, 
and the fluid layer in between, see Fig. 1. In the case that the two parallel walls have identical 
temperatures, the dissipative work of the viscous force will still lead to a parabolic temperature 
distribution inside the fluid layer. From the viewpoint of dimensional analysis, the heat transfer 
of such a system shall be governed by two dimensionless parameters, namely the Prandtl number 
and the Eckert number. As the Prandtl number of the lattice BGK modeled fluid has an invariable 
unit value[12], the Eckert number, defined as 

E=— ^ -, (13) 

becomes a decisive parameter. The definition of the Eckert number tells that it is a measure of 
ratio between the heat due to friction and the heat conducted by temperature difference. Note 
that this number will become infinite when the temperatures of two parallel walls are the same. 
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Figure 1: Velocity and temperature distribution in the Couette fiow. 

The analytical solution of the steady temperature distribution in the transverse direction of 
the channel can be obtained by directly solving the Navier-Stokes equations, which read as. 



Ti=To 
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In the following calculations, numerical results will be normalized by (Ti — Tq) if Ti ^ Tq or by 

otherwise. The linear distribution of the fiow velocity will not be checked here, as it can 
always be obtained as long as the steady stage of fiow is reached. 

3.2 Results of Qian's model(2i)13yg) 

The model used here is a typical lower order thermal lattice BGK model on the composed 
square lattice. Conditions for the calculation are set as follows: lattice size 64 x 32, time step 
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Figure 2: Temperature profiles in tlie Couette flow under different Eckert numbers. Solid lines are 
analytical solutions and symbols of crosses are the numerical results calculated by using Qian's 
thermal lattice BGK model. 



10,000, moving boundary at j/ = 1 and non-slip and fixed temperature boundaries. The analytical 
solutions and the numerical results are shown in Fig. 2, from which one may conclude that the 
numerical calculation is accurate only when the Eckert number is small enough, that is, when 
the temperature distribution is mainly controlled by the heat conduction between the two walls. 
Notice that the Mach number Ma = where a,, is the adiabatic sound speed, was kept to be 
smaller than 0.05 for all the three cases, in order that deviations caused by higher order terms 
contribute diminishing effects. It is clear, because of this low Mach number condition, that errors 
occurred under the large Eckert number are solely due to the lack of an exact expression for the 
viscous work in the energy equation. 

It is crucial to derive the concrete expressions of the hidden terms in the energy equation of the 
modeled fiuid, since one may then know the source of those large errors and improve the quality of 
lower order modeling if possible. We have derived the structures of these terms under two optional 
constraints. 



pk pk 

^pkjpkO = , ^ ^pkjpkl = . 

pk pk 
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As mentioned before, these optional constraints are used to minimize the anisotropic effects. The 
hidden terms can be explicitly written out in five parts[12]. 



Parti : -da{pu'^Ua) 



Part2 : da[{£,l + £,26aj3)dj3{peUaUj3) + {£,3 + £,46al3)dj3{pUaUj3)] 
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Here, we use A^'*'^] to represent the sum of permuting products of the Kronecker tensors, namely 
(6ai3^-y6(^ + • • •) appearing in Eq. (4). The nonlinear response coefficients are defined in the 
foUowine formulas, 



(19) 
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As the definitions of (i and (2 involve ^pk, the contribution of PartA is anisotropic. Observing 
other parts, we see that there always exist terms of du^ order except in Parti and Part5, which 
resulted from the higher order terms in the Euler level energy equation and the Navier-Stokes level 
momentum equation. Since the du"^ terms have the same order as the normal viscous work terms, 
their combined infiuence should be the direct source of errors under the large Eckert number. Such 
errors, as discussed so far, are independent of the Mach number. 



3.3 Results of a revised lower order model(2i)13yC) 

The quality of the lower order thermal lattice BGK model can be improved by making changes 
in the optional constraints stated in the previous section. Actually we may relax one of such 
constraints, Eq. (16), and impose two new constraints as follows 
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(20) 
(21) 



With the use of the definition of the shear viscosity in Eq. (11), it can be shown that these 
constraints will lead to the cancellation of terms of du"^ order appearing in Part2 and Parti of 
Eqs. (18). Nevertheless the relaxation of Eq. (16) will bring about an additional part of anisotropic 
errors, which involves the lattice symmetric parameter Ap^ and the higher order Kronecker tensor 
'^afiiiCi- Note that the anisotropic errors did not occur in the previous calculation, because 
ilpfcs' are zero in the two dimensional space. However, Apj,s' are non-zero parameters in all the 
dimensions, so that anisotropic errors will play their roles in the numerical calculation using the 
revised lower order model. 
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Figure 3: Change of the underlying lattice in the revised lower order model. The 2D13VQ model 
uses 00, 11, 12 and 21 sub-lattices and the 2D13VC model uses 00, 12, 21 and 22 sub-lattices. 



The satifaction of Eq. (21), may not be realized without a further consideration of the lattice 
geometry. Since 00, 11, 12 and 21 sub-lattices are employed in the 2D13VQ model, 02i becomes 
the only non-zero parameter in the 6th rank velocity-moment tensor, and so is the <^2i in the 4th 
rank tensor. Now that j2io and j2ii would have already been decided by the cppks' constraints, 
there would be no room for them to satisfy those 0pfcs' constraints. The solution is to replace the 
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Figure 4: Temperature profiles in the Couette fiow obtained by using the revised lower order 
thermal lattice BGK models. Solid lines are analytical solutions, others are the results of numerical 
calculations. 

11 or 12 sub-lattice with the 22 sub-lattice. Parameters such as 02i and 622, f2\ and (^22 all have 
non-zero values, so that the satisfaction of both the <,Cpfcs' and 0pfcs' constraints becomes possible. 
We employed 00, 12, 21, 22 sub-lattices for the 2D13VC model, as shown in Fig. 3, taking care of 
the numerical stability. If the 11 sub-lattice is used instead, the difference of particle populations 
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on different sub-iattice wiii be so iarge that particies may get neariy depieted on some finks of 
certain sub-fattice. In that case, even a very smaff fluctuation could cause the appearance of 
negative particle distributions and thus the numerical instability is triggered on in the calculation. 

Calculations were performed under the conditions described before. The results of the 2D13VC 
model and the analytical solutions are compared with each other in Fig. 4. It is found that errors 
have been greatly reduced even when the Eckert number becomes infinite. The remaining errors 
can be recognized as two parts. The first part of errors is brought by terms of or higher order. 
These errors become obvious when the Mach number is too large, see the broken line in Fig. 4. 
The loss of the symmetry for this line is due to the different local Mach numbers distributed across 
the section. Another part of errors comes from the anisotropic terms. The behavior of these errors 
in this case is an oscillating one, but small in both scale and magnitude. We are aware that even 
the anisotropic errors can be completely deleted if both the 11 and 12 sub-lattices are employed 
in the 2D13VC model. But this could make the lower order model cost the same amount of 
the computer resource as the higher order model, which is absolutely unfavorable, especially in 
the three dimensional case where the lower order model could save nearly half of the computer 
memory consumed by the higher order one. 

3.4 Results of the higher order model(2i)16y) 

It is clear now that the lower order thermal lattice BGK models suffer from different errors under 
different conditions. Although they are improvable, errors cannot be completely avoided. On the 
other hand, the higher order model was designed to eliminate all these deviations. It was proved 
before[13] that the deviation terms existing in the macroscopic momentum equation was removed 
for the higher order model . The calculation carried out here further shows that the modeling of 
heat transfer by using such a model is also advantageous. From the comparison of the numerical 




Figure 5: Temperature profiles in the Couette fiow obtained by using the higher order lattice 
BGK model. The solid line is analytical solution and various symbols are numerical results under 
different Mach numbers. 

results and the analytical solutions shown in Fig. 5, we may know that the higher order model 
behaves extremely well no matter how the Eckert or the Mach numbers change. 
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4 Conclusion remarks 



Through theoretical and numerical analyses presented above, the ability of the thermal lattice 
BGK models in modeling heat transfer in fluid flows was well demonstrated. Furthermore, we 
may contribute some advice to the users of such models: If one wants to save computer resource, 
specially in those 3D situations, and does not care about small scale anisotropic errors, it is safe 
to use the lower order thermal lattice BGK models in the incompressible regime. On the other 
hand, if high accuracies are required, then the higher order model is preferable. When flows enter 
the transonic regime, it is essential to use the higher order model to avoid large deviations in the 
simulation results. 
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